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As a step toward analyzing second-harmonic generation (SHG) from crystalline Si nanospheres in glass, we 
develop an anisotropic bond model (ABM) that expresses SHG in terms of physically meaningful parameters 
and provides a detailed understanding of the basic physics of SHG on the atomic scale. Nonlinear-optical (NLO) 
responses are calculated classically via the four fundamental steps of optics: evaluate the local field at a given 
bond site, solve the force equation for the acceleration of the charge, calculate the resulting radiation, then 
superpose the radiation from all charges. Because the emerging NLO signals are orders of magnitude weaker 
and occur at wavelengths different from that of the pump beam, these steps are independent. Paradoxically, the 
treatment of NLO is therefore simpler than that of linear optics (LO), where these calculations must be done 
self-consistently. The ABM goes beyond previous bond models by including the complete set of underlying 
contributions: retardation (RD), spatial-dispersion (SD), and magnetic (MG) effects, in addition to the anhar- 
monic restoring force acting on the bond charge. Transverse as well as longitudinal motion is also considered. 
We apply the ABM to obtain analytic expressions for SHG from amorphous materials under Gaussian-beam ex- 
citation. These materials represent an interesting test case not only because they are ubiquitous but also because 
the anharmonic-force contribution that dominates the SHG response of crystalline materials and ordered inter- 
faces vanishes by symmetry. The remaining contributions, and hence the SHG signals, are functions entirely of 
the LO response and beam geometry, so the only new information available is the anisotropy of the LO response 
at the bond level. The RD, SD, and MG contributions are all of the same order of magnitude, so none can be 
ignored. Diffraction is important not only in determining the pattern of the emerging beam but also the phases 
and amplitudes of the different terms. The plane-wave expansion that gives rise to electric quadrupole/magnetic 
dipole effects in LO appears as RD here. Using the paraxial-ray approximation, we reduce the results to the 
isotropic case in two limits, that where the linear restoring force dominates (glasses), and that where it is absent 
(metals). Both forward- and backscattering geometries are discussed. Estimated signal strengths and conversion 
efficiencies for fused silica appear to be in general agreement with data, where available. Predictions are made 
that allow additional critical tests of these results. 



I. INTRODUCTION 

Second-harmonic generation (SHG) is becoming an in- 
creasingly important diagnostic tool for a wide range of ap- 
plications. It is a particularly important probe for studying 
planar interfaces between centrosymmetric crystals and over- 
layers with randomly directed bonds, because it is dipole- 
allowed only at the interface where the bonds are simultane- 
ously asymmetric and well ordered. 

Recently, Figliozzi et alii found that SHG signals gener- 
ated in transmission from crystalline Si nanospheres (nSi) dis- 
persed in glass were enhanced significantly when driven by 
two beams with crossed polarizations. Enhancement of any 
nonlinear-optical (NLO) signal is automatically of interest, 
because in principle NLO signals contain significantly more 
information about materials systems than the linear-optical 
(LO) response, yet are intrinsically much weaker. SHG from 
the dispersed-nSi configuration was recently analyzed from 
the macroscopic perspective by Brudny et al. 2 and Mochan 
et al.^ in the former case for a single isolated nanosphere 
and the latter for arrays of nanospheres. These authors used 
the "dipolium" approximation^ where the inclusions and host 
are described macroscopically by linear, isotropic dielectric 
functions. The far-field SHG response was obtained by cal- 
culating the effective dipole of the inclusions as a spherical- 
harmonic expansion of the internal and external fields of 
a given nanosphere, then applying standard radiation equa- 
tions. Various observations were explained, for example the 
(E ■ V)E symmetry of the SHG intensity, its dependence on 



sphere size, the importance of screening in determining the 
contributions from the interiors of the nanospheres, the emis- 
sion of SHG radiation in a cone for disordered dispersions 
of nanospheres, and the relatively small intensity of the SHG 
signal from glass. 

While macroscopic treatments efficiently distinguish be- 
tween allowed and forbidden contributions, they are unable 
to relate allowed responses to atomic-scale parameters, or to 
provide the same level of understanding of the different con- 
tributing processes. In particular, the following questions still 
need to be answered: (1) how does the SHG intensity from the 
nSi inclusions compare to that from planar Si-Si02 interfaces; 
(2) what are the relevant parameters; (3) what is the maximum 
intensity that can be obtained; and (4) is this maximum signal 
useful, or simply given by a combination of already known pa- 
rameters? While much larger SHG signals might be expected 
from dispersions of nSi inclusions in a transmission configu- 
ration simply because the interface area greatly exceeds that of 
a planar interface, the larger area is offset by the fact that the 
first-order anharmonic SHG signals from the opposite sides of 
the nanospheres cancel. Therefore, the signal is proportional 
to the gradient of the driving field instead of the field itself. 2 
In addition, contributions are limited in depth to the coher- 
ence length in the material. Finally, there is the question of 
whether the enhanced SHG signals observed with dual-beam 
excitation provide useful information. The atomic-scale mod- 
eling done below shows that the contributions of the three 
underlying mechanisms, retardation (RD), spatial dispersion 
(SD), and magnetic (MG), can all be predicted from the LO 
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response and beam characteristics, hence do not necessarily 
provide new information even though improved geometries 
may generate large signals. 

In addressing these issues we found it necessary to ex- 
tend our previous simplified bond-hyperpolarizability model 
(SBHM), which expands on the even simpler isotropic force 
model discussed for example in Shen. 5 In the SBHM, SHG 
is expressed as radiation arising from the anharmonic mo- 
tion of charge localized in bonds assuming that the only mo- 
tion relevant to SHG is that along the bond direction itself. 
The SBHM successfully describes, with many fewer parame- 
ters than previously required, a wide range of NLO phenom- 
ena including SHG& and FHGi from Si-insulator interfaces, 
dipole-forbidden SHG^ and THG^ from centrosymmetric ma- 
terials, and the generation of THz radiation from III-V semi- 
conductor surfaces*^ In addition, the parameters are physi- 
cally meaningful, and by incorporating crystal symmetry at 
the atomic level, macroscopic tensor properties are obtained 
automatically. However, as recently, and correctly, noted by 
McGilp^ the SBHM has limitations regarding quantitative 
interpretation. Given the simplicity of the approach this is not 
surprising, but it needs to be explored further. This is a second 
objective of this work. 

Accordingly, in the present paper we generalize the SBHM 
to a more complete description, the anisotropic bond model 
(ABM), which includes charge motion transverse to the bond, 
RD, SD, and MG effects, including SD effects arising from 
beam geometry, and SHG signals for off-axis observation an- 
gles, i.e., the role of diffraction. In developing our expressions 
we follow the approach of Peng et al.^ framing the calcu- 
lations in terms of the fundamental 4-step process of optics: 
(1) evaluate the local field at any given charge site that re- 
sults from the driving (source) field; (2) solve the mechanical 
equation F = ma to obtain the acceleration of the charge; (3) 
calculate the radiation that results from the acceleration; and 
(4) superpose the radiation from all contributing charges. For 
random media we show that step (4) factors into two parts: 
(4a) average over all possible bond orientations at a single 
site, then (4b) calculate the properties of the emerging beam 
by Fourier-transforming the envelope function of the incident 
radiation. Although not required here, if appreciable energy 
were transferred from the driving to the generated beams, then 
it would be necessary to (5) evaluate the energy extracted from 
the pump beam as a function of position, with the subsequent 
correction of the local fields evaluated in step (1). We find 
that for random materials the RD, SD, and MG contributions 
are all of the same order of magnitude and must all be con- 
sidered. Finally, all aspects, including off-axis observation 
and diffraction, combine to yield a much richer SHG response 
than previously assumed. 

Aside from including bond anisotropy and additional mech- 
anisms, the approach is essentially the NLO equivalent of that 
which Ewald 12 and Oseeni^ used nearly a century ago to de- 
rive the Ewald-Oseen theorem of LO. Paradoxically, from this 
perspective NLO is simpler than LO. In LO the radiated fields 
have the same wavelength as the driving field and similar in- 
tensities, so steps (1), (3), and (4) must be evaluated self- 
consistently. In contrast, for NLO the radiated fields are typ- 



ically orders of magnitude weaker than the driving field and 
occur at different wavelengths, so all steps are effectively in- 
dependent. This allows NLO calculations to be done sequen- 
tially, to levels of approximation that are also independent and 
may be adjusted to meet particular requirements. This is one 
of the few cases where a nonlinear problem is simpler than its 
linear equivalent. 

Advantages of an atomic-scale formulation include a bet- 
ter understanding of the physics involved. In this classical 
model NLO is a result of distortions of the nominally sinu- 
soidal waveform of the emitted radiation reaching the ob- 
server. The obvious contributing factor is anharmonic motion 
of a charge. This can be due to an anharmonic restoring force 
(intrinsic anharmonicity), spatial nonuniformity of the driv- 
ing field (spatial dispersion), or the magnetic field associated 
with the driving wave. With respect to acceleration, there can 
clearly be no distinction between anharmonic motion result- 
ing from an anharmonic restoring force, a field that is slightly 
larger at one limit of the excursion than the other, or a force 
that is velocity-dependent. All these effects enter in step (2). 
However, another source of distortion is the finite speed of 
light. This causes the signal reaching the observer from the 
far limit of the excursion to be delayed slightly relative to 
that from the near limit, resulting in a waveform distortion 
equivalent to phase- or frequency-modulation. The retarda- 
tion contribution enters in step (3). Retardation involves the 
same first-order expansion of a plane-wave factor that leads to 
the electric quadrupole/magnetic dipole contribution of LO, 
but the physics is quite different. This mathematical similarity 
has led to confusion in the past, and we clarify the distinction 
below. 

Taking into account the complete set of mechanisms con- 
tributing to SHG became a larger project than expected, so in 
the present paper we restrict applications to single-beam ex- 
citation of disordered materials and reduction of the resulting 
expressions to the case where the bond charges are isotrop- 
ically polarizable in LO. We discuss two limits, first where 
the restoring force dominates the acceleration term (glasses), 
and second where the restoring force is absent (metals). This 
reduction, done in the paraxial-ray approximation, highlights 
the roles of the different underlying mechanisms, the differ- 
ence between forward- and backscattering configurations, and 
allows a simple expression for signal strength and conver- 
sion efficiency to be obtained. The present work represents 
a necessary first step toward our goal of understanding, at 
the atomic level, SHG from Si nanoinclusions in glass under 
crossed-beam excitation, and is further justified by the fact 
that disordered materials are ubiquitous in many fields. 



II. AMORPHOUS MATERIALS 
A. Fields at bond sites 

In this section we consider step (1), define basic quanti- 
ties, and discuss the connection between first- and second- 
harmonic fields. We suppose that the relevant quantities are 
electrons of charge q = — e located in bonds j at positions 
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r q = Fj + Arj , where the fj are the equilibrium positions of 
the charges relative to the origin of a coordinate system de- 
fined in the laboratory and the Ar) are the displacements that 
result from the time-dependent forces acting on them. We rep- 
resent the directions of the bonds by bj, where for Si-O bonds 
the bj point from Si to O. 

For amorphous materials that are homogeneous on meso 
and macroscopic length scales, the driving field can be as- 
sumed to be approximately a plane wave with frequency ui, 
envelope function E (f q ), and wave vector k Q = {um\/c)k , 
where ri\ is the refractive index of the material at w. We as- 
sume that the Fresnel reflectance coefficients have been ap- 
propriately taken into account at the surface of the material to 
yield the correct amplitude E a of E (f q ) within the medium. 
Then the field at the j th charge can be written to first order in 
Afj as 



E(r q ,t) = E (f q )e ik ->^- iut 



E Q {r 3 + Af J )e i 'fi +A? ^- wt 
[l + Afj-V^Eoir^e^-^ 
[1 + Afj • V^]4e -ihrt . 



(1) 



For clarity in the following equations, we let Ej = 

E (fj)e lk °' Vj contain the spatial dependences of the envelope 
and phase. The SHG nature of the correction term follows be- 
cause Arj is also proportional to Ej, as shown below, so the 
gradient term nominally has a time dependence e~ l2uJt . 

However, the coefficient of a e - * 2 "* term is not simply the 
product of the coefficients of the parent e~ luJt terms, but must 
be reduced by a factor of 2 for the following reason. Ob- 
servables are real quantities, so e~ luJt = cos(uit) + i s'm(ojt) 
is actually shorthand for Re(e~ lut ) = cos(arf). Thus the 
product of two e~ luJt terms is really a product cos 2 (wi), 
sin(wi) cos(wi), or sm 2 (cjt), or some combination depend- 
ing on the phases of the parent coefficients. All trigonometric 
identities taking ut products into 2u>t forms involve a factor 
of (1/2). We introduce this factor in the far-field radiation ex- 
pression Eq. (|26V We retain the e~ l2uJt notation so average 
intensities can be calculated in the usual way. 



B. Force equation 

The general form of the force equation for SHG is 

d 2 Af(t) 



F = ma = m- 



dt 2 



qE(f,t)+q- x B(f,t) 
c 

- k 1 ■ Af(t) - k 2 ■ -Af(t)Af(t), 



(2) 



where m is the mass of q and ki and k 2 are second- and third- 
rank tensors describing the linear (Hooke's Law) and first- 
order anharmonic restoring forces, respectively, v = dArJ dt, 
and the magnetic-flux density B(f q ,t) associated with the 



driving field is B(r q ,t) = —(ic/u))Wp x E(r q ,t). In con- 
trast to the SBHM, we do not assume the force equation to be 
one-dimensional. To find the displacements Arj, we substi- 
tute Eq. (Q]i into the force equation to obtain 



in 



d 2 ^m 
dt 2 



= q[l + Ar j (t)-V? j ]E j e- i " t 

q dArAt) : ,, 

c dt 

- hi ■ Arj{t) - k 2 ■ ■Af ? i (f)Ar*,-(i), 



(3) 



where Bj = (— ic/w)V 'p x Ep j . From the form of E(f q , t), 
we can assume that 

Afj(t) = Afye"^* + Ar 2] e- l2 <»\ (4) 

where Ar\j and Ar 2 j are time independent. Substituting this 
expression in Eq. (01 yields 

- muu 2 Ar l3 e~ iult - Amuj 2 : Ar 2j e~ i2uJt 
= q(l + (Ar lj e- i " t )-Vp j )E j e- i " t 

- g(Arij-e- iwt ) x (V rj x Eje~ iut ) (5) 

- hi ■ (Af 1] e-^ t + Ar 2] e- a ^) 

- k 2 ■ •Af y Afi i e- i2wi . 

Since Afj is at least first-order in E, the magnetic term is at 
least second-order in E. Since we are only concerned with 
SHG, we neglect terms of order (3w) and (4cj), which would 
contribute to THG 9 and FHG^ respectively. 
Isolating the first-harmonic terms we have 



mLO Afij = qEj — k\ ■ Arij. 



(6) 



While this can be solved in general by matrix methods, we 
now introduce the approximation that k\ and k 2 are diagonal 
in the local coordinate system of the bond, where the z axis is 
defined by the unit vector b parallel to the bond. Diagonaliza- 
tion is equivalent to assuming that the bonds are rotationally 
symmetric. While bonds in some systems are not rotationally 
symmetric, we make this simplifying assumption to elucidate 
the underlying physics. Obviously, if desired all tensor com- 
ponents of the restoring forces could be kept. 

We also define the unit vector i, which is perpendicular to b 
and lies in the b — E plane. Thus t is given by 



t 



E-b(b-E) \ A E 2 - [b-E 



Then ki and k 2 can be written as 

ki = bbnu + Unit, 

k 2 = bbbK2i, 



(7) 



(8) 



(9) 



where Ku and k 2 i are the longitudinal linear and anharmonic 
restoring-force coefficients, respectively, and Kit is that for 
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transverse displacements. With the assumption of rotational 
symmetry, K 2t does not exist. However, transverse contribu- 
tions are still possible through the RD, SD, and MG terms. 

Substituting these expressions into Eq. (0 and taking dot 
products with b and t leads to the two first-order equations 



Af w = Ar^bj = q(bj ' Ej \ b r , 



Af ljt = Arijtij 



,2 J 



U. 



(10) 



(11) 



(12) 



Kit — rnoj 

Repeating the process for the second-order terms leads to 

Ar 2 ji = Ar^jibj 

= q{Ar ir W ?j )(bj- 4) 

- g(Afij x (V r - x E 3 )) ■ bj 

- K 2 iArijiAriji bj/(nu - 4mw 2 ); 

Ar 2jt = Ar 2 jdj 

= q(Af f lj ■ V fi ){f s -Ej) 

- g(Afy x (V^. x E 3 )) ■ ij ij/(K U - Amiv 2 ); 



Af 2 j (bxt) = Ar 2 j{txb){bj x ij) 



(13) 



(14) 



(Anj x (V r - x Ej)) ■ (bj x tj) 
x (bj x ij)/(nit — 4mw 2 ). 

Equation ( TBl i is necessary because the magnetic force gener- 
ates a component that is perpendicular to both b and t. Equa- 
tion (f]~2b shows that there is no qualitative distinction between 
the intrinsic anharmonicity of a bond and an anharmonicity 
generated by a field, as expected. These are the expressions 
from which the acceleration, and therefore the far-field signal, 
will be calculated. 



C. Far-field radiation from accelerated charges 

We now consider step (3). We follow the development of 
Peng et al., 8 but take into account explicitly the reduction in 
propagation speed caused by refractive indices n v that are 
different from 1. The two that need to be considered are 
ti(oj) = \J t(uj) = n\ for the incoming wave and n(2uj) = 
y 7 e(2ui) = n 2 for the emitted SHG radiation. Accordingly, 
we write k a = k k = (umi/c)k and k = kk = (2um 2 /c)k 
for the incident and emerging radiation, respectively, where k 
points in the direction of the observer. 



The general expression for the four-potential of an acceler- 
ated point charge in the medium in Fourier-component form 
is 



cf>(r,t),A(f,t) 



d S r'dt' 



cp(f',t'),J(r',t') G v (r,r\t,t') 



(15) 



where 



„ 5(t-t' 
G„(r,r',t,t') = ^ 



"0 



\f — r'\ 



(16) 



is the Green function, p and J = pv are the charge and cur- 
rent densities, respectively, and n v is ri\ or n 2 according 
to whether the frequency of interest is u> or 2u>. Here, p 
and J are associated with the j th point charge q located at 



Afj(t). Then 
P j(r',t')=q5(f'-rj-A?j(t')): 



, ^ dAfj(t') 



Jj(f',t')= Pj (rV) 



dAfj(f) 
dt' 



dt' 



6(r' -?j-Afj(t')). 



(17) 



(18) 



The far-field (r, t) that results from q is given by 



c dt V ^,t) 

l dl jx (f,t) 

c dt 



(19) 



where Aj±(r, i) is the component of Aj(r, t) that is perpen- 
dicular to the line between the origin rj + Arj of the radia- 
tion and the observer at r. The second line of Eq. (T% follows 
because V(j>j in the first line removes the longitudinal com- 
ponent of Aj, leaving a purely transverse potential. Thus we 
need evaluate only Aj(f,t). This can be obtained relatively 
simply because Aj is already of first order in v /c, where v is 
the velocity of q. 

Substituting Eq. ( fT8l into Eq. ( fT3T > and performing the inte- 
gration over r 1 yields 



A,(f,t) 



q f ^ dArj (?) S(t-t'-^\r- r, Afj (t')\) 



dt' 



\r - r 3 - Ar 3 (t')\ 



(20) 



The integration over t 1 is nontrivial because Ar) (t ') is also a 
function of t' . However, to first order in 1/c we can expand 



^ |f- rj - Arj(t')\ a M - ■ fj - 'fk ■ Afj(t>), 

(21) 



c 
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5[t-t'- — \r-fj - Afj(t')\ 



(22) 



the third is a combination arising from the spatial dependence 
of the field (SD, MG) and the intrinsic anharmonicity of the 
bond (the third term in Eq. (fT2l). Because the RD contribu- 
tion originates in propagation, not acceleration, the use of the 
common expression 



wherei = t—n v r/c+n v k-fj/c. This is still a self-consistent 
expression, but to first order in 1 jc we can substitute t a for t' 
in the argument of Afj (t 1 ). We obtain finally 



t' = t ret ~to + —k ■ Afj (t a ) , 



(23) 



where t ret is the retarded time. The integral over t 1 can now 
be performed, and we obtain 



(24) 



Substituting Eq. into Eq. (124-b yields the contribution from 
the j th charge: 



Aj(f,t) = (Afye"^' + 2Af 2l e-^ t ' 



rc 
iujq 
rc 



Afi f e _ ***o-Afij e -ik k-rj e ik a r-iujt 



+2Af 2j e- ik ^e ikr - l2uJt 



- — Afi a e~ ikk ° *i e ikoT ~ iut 
rc 



UJ 2 qn 2 



r Af lj (k ■ Af l3 )e- %k ^ e ikr ~ i2u]t 



rc 

i -^-Af 2 Ae- ik ^e ikr - i2ut 
rc 



The far field signal Ey then follows from Eq. (TT9b : 



(25) 



E ff (f,t) 



I -kk 



2 

lo q 



Ar\j e 



— ik k-fj gik r — iut 



3 __ 2 Afi-, (fc • Ar*ij)e 



— ik-fj ^ikr—i2tut 



ik-rj gikr—z2ujt 



(26) 



Here, / — kk is the projection operator that eliminates the 
longitudinal component and hence performs the function of 
—V(j>j. As with — V<f>j, I — kk does not affect the orthogonal 
component, which will be found to be significant when we 
discuss term cancellations in Sees. Ill B and III C. In the two 
nonlinear terms of Eq. (|26| |. we have now incorporated the 
factor of (1/2) associated with the change of time dependence 
from (e~ luJt ) 2 to e ~ l2uJt as discussed in Sec. Ill Al 

Equation d26l i is a general expression for linear and second- 
order far-field radiation from a moving charge in terms of dis- 
placements from its equilibrium position. The first term is the 
linear response. The second term is the RD contribution, and 



Eff 



1 d 2 a ± 



(27) 



leads for this term to an error of a factor of 2. 



To address a point that has caused difficulty in the 
past, we note that the RD term above and the electric 
quadrupole/magnetic dipole (EQ/MD) terms of LO both re- 
sult from an expansion of a phase term e lk ' r to first order in 
k ■ f. However, the physics, and consequently the nature of 
EH , is different in the two situations. In LO p(r,t) is as- 
sumed to be a moderately extended but stationary charge den- 
sity with a multiplicative time dependence e~' luJt , thus hav- 
ing the form p(r, t) = p (f)e~ %u}t . Here, the t' integration is 
trivial but the f ' integration is not. Performing the t! integra- 
tion yields a multiplicative time factor e~ lut and a phase term 
e r that is part of the electrostatic Green function. Because 
the current J needed to calculate A has no obvious represen- 
tation in this case, appropriate vector-calculus identities are 
used to convert the integration of J into a first-moments in- 
tegration of f'p(f')w^ The dipole approximation follows by 
taking e = 1, with higher-multipole moments generated 
from higher-order expansion terms Thus the LO expan- 
sion gives rise to multipole moments but no higher harmonics. 

In contrast, in the present work p(f, t) describes a moving 
point charge qS(f ' — r (t)), where f (t) = fj + Afe~ lut . As 
seen above, the f ' integration is now trivial but the t' integra- 
tion is not. We obtain here higher harmonics but no multipole 
moments. Thus what Peng et al. labeled EQ/MD in refi£ is 
due to retardation. That in refi^ is actually due to spatial dis- 
persion. 



D. Superposition of radiation; averaging and diffraction 

In the following we assume that the charges are driven co- 
herently, so fields must be added rather than intensities. This 
is expected, and the validity of the assumption demonstrated 
experimentally by the vanishing of SHG for amorphous mate- 
rials in the forward direction. 

Returning to Eqs. (Q}, (fT0b-(fT3]l, and d26b , the fj depen- 
dence of the SHG signal is either E 2 (fj)e^ 2ko ~ k ^' rj or a 

derivative of the form E a (fj ) (dE a (fj ) / dx)e l( ~ 2k ° ~ k ^ p i . Both 
are slowly varying on the atomic scale, whereas the bond di- 
rections b and t vary essentially randomly from site to site. 
Given this large difference of scale we can factor step (4) into 
two parts: averaging over bond orientations, effectively at a 
single site; then evaluating the sum over all fj. 
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1. Bond averages 

We consider first averaging over bond directions. This is 
accomplished by writing 



6 = b x x + b y y + b z z 
= xsinOcoscj) + ysinOsincj) + zcos6 



then performing the operation 



(28) 



(29) 



The calculation is simplified by grouping the products in- 
volving bond directions into dyadics, triadics, etc., then con- 
sidering symmetry. For example for LO the bond averages 
that need to be evaluated occur as dyadics 66 and it. In the 
Cartesian-coordinate representation 66 has 9 terms b x b x xx, 
b x byxy, etc., but only 3 survive the averaging process because 
any component involving an odd number of projections aver- 
ages to zero. 

By this reasoning the triadic 666 clearly vanishes identi- 
cally, so by Eq. (TlZt there can be no k^i contribution to SHG 
in amorphous materials. Not surprisingly, both microscopic 
and macroscopic considerations therefore lead to the same 
conclusion. However, the implications here go further. The 
terms that remain are functions only of the LO response and 
the configuration geometry, so the amount of new information 
obtainable by SHG in amorphous materials is limited to the 
separation of longitudinal and transverse components of the 
LO response, no matter what geometries are used to enhance 
the SHG signal. 

We now consider the RD contribution to SHG. The terms 
that need to be considered are 6666, bbtt and permutations, and 
tttt. If desired, the results can be decomposed into irreducible 
tensor representations, although we do not do this here. In the 
calculations that follow we take advantage of the absence of 
a preferred direction in amorphous materials. Hence without 
loss of generality we assume that k a = k a z and E a = E Q x. 
The expression to be evaluated is then 



I-kk 



■ — [ dtt (Aru ■ k ' 
4W V J (30) 



x Af lje ^ k ^e lkr ^ l2uJt . 



The result is 



ff _ .LQ qn 2 



RD = ~i 



15rc 3 



I-kk 



x(k • i)(3Cf + ACiCt + SC 2 ) 

+ (y(k ■ y) + z(k ■ z)) (Q - C t ) 2 



i(2k —k)-rj 



Akv—ilujt 



(31) 



where to simplify the expression we define 

9^9 



a = 



a 



Ku — raw' 



Kit — muj' 



(32) 



We write the x component etc. of k as k(k ■ x)x so we can 
move the magnitude of k to the prefactor and therefore elimi- 
nate an easily overlooked source of error. The separate longi- 
tudinal and transverse contributions can be obtained by setting 
Ct = or Ci = 0, respectively. The sum over fj is clearly a 
Fourier transform of the square of the envelope function of the 
driving beam, and will be evaluated in Sec. lIID2l 

We consider next the contributions from SD. These are 
given by 



E 



ff 

SD,j 



2oj 2 q 



I-kk 



4tt 



q(Ar ll -V rl )(b-E J ) i 
ku — Amu 2 



(33) 



Kit — 4mw 2 



We evaluate Eq. (1331 by dividing the field gradient into lon- 
gitudinal and transverse parts with respect to k , i.e., letting 

V^E(fj,t) = x^ + y^j + izk Q . After performing the 
averages we obtain 

->ft , 2uj 2 q 



15rc 2 



x < x 



.dEjfj) 
dx 



7;A + 2C,A + 2CtA + 8C t D t ) 
(d - C t ) (D t - D t ) 



dy 



^ gi(2k — k)-rjgikr—i2u>t 



where 



D, = 



Ku — Amu 2 ' 



Ku — Amu; 2 



(34) 
(35) 



The envelope function for the z component is the same as that 
for the RD contribution, but those for x and y involve gradi- 
ents of the driving field. 

The MG contribution is given by 



E 



if 

MG,j 
X 



2uj 2 q 



I-kk 



1 

4tt 



dil {A UAfij x (V Fj xEj))-b 
+ D t [(Afy x (V r - xEjj)-t 



t 



+A 



(Afy x (Vfj x Ej)) ■ (6 x i) (6 x £)| 



^, ^ — ik-rj^ikr — i2ujt 



(36) 
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Here, all three dimensions are involved. By suitable vector- 
calculus identities the double-cross-product operation can be 
cast into apparent spatial-dispersion form, (E ■ V)Er but an 
exact cancellation of the resulting dominant terms makes this 
approach unproductive. After performing the cross-product 
operations with the assumed propagation and field directions 
and then averaging over bond orientations we obtain 



2u?q 



3rc 2 

x (2C t A + QDt) 



I-kk\ -J2 E °( r 
dE {r,) , 



dy 



izk E {fj) 



x e 



i(2k —k)-rj ikr—i2ujt 



(37) 



pie the result is 



nNk x W 2 E 2 
' A{2k -k z ) 



g \ x ' y I 



(41) 



This is also a Gaussian beam, but with a nodal line passing 
through the center. This is the analytical representation of the 
two-lobed pattern reported by Figliozzi et alii for various con- 
figurations of SHG from amorphous material and spherical Si 
nanoinclusions. 



E. Net results 



2. Diffraction 

We consider now the sums over fj . These not only yield the 
geometric properties of the emerging SHG beam, but also af- 
fect the phases and amplitudes of the prefactors of the individ- 
ual constituents. In the derivation below we assume forward 
scattering, but will discuss backscattering in Sec. HHP I We 
consider throughout only single-beam excitation. Crossed- 
beam configurations follow the same principles but are com- 
plicated by the need to consider large observation angles, so 
will be treated in a subsequent paper. 

We assume that the incident beam is Gaussian. For the RD, 
z SD, and z MG contributions, the relevant sum is 



E 



E 2 e -2(x 2 +y 2 )/W 2 



:> i(2k — k)-rj 



(38) 



where W is the width of the incident beam and for our con- 
figuration k Q = k a z. Converting the sum to an integral we 
have 



-» N [ dxdy ( dz, 

J-oo JQ 



(39) 



where A*" is the volume density and L the thickness of the 
sample. The integrals are all standard and we find 



nNW 2 E 2 
2(2Aj k z 



(40) 

where we have assumed that L is much larger than the coher- 
ence length l/(2fc — k z ). As expected, the emerging beam 
also has a Gaussian cross section, with a contributing volume 
determined by the size of the original beam and the coherence 
length of the configuration. 

For the x and y SD components and the y MG component 
the integrals are also standard. Taking the x term as an exam- 



We now combine the results of the above sections. The 
overall RD contribution is 



E 



ff 



iruj 3 qNW 2 E 2 



I-kk 



30rc 3 (2k o -k z ) 

x(k ■ x)n 2 (3Cf + ACiCt + 8Cf) 

+ (y(k ■ y)n 2 + z{k ■ z)n 2 J (Q - C t f 



x e 



-(kl+k 2 )W 2 /Sikr-i2uit 



(42) 



The corresponding expressions for SD and MG are respec- 
tively 



E"(f,t) 



noj qNW E 2 
15rc 3 (2k - k z 



I-kk 



x(k ■ x)n 2 (3CjDj + 2QA + 2C t D l + 8C t D t 
+ (y(k ■ y)n 2 + zn x ) (Q - C t ) (D l - D t ) 

x e 



(43) 



E 



ff 



ttlo qNW E 2 



2p2 r 



I-kk 



(44) 



3rc 3 (2fc -fc 2 , 
y(k ■ y)n 2 + zni] (C t D t + 2Cj D t ) 

x e -(kl+k 2 y )W 2 )/B e ikr-i2ujt _ 

Equations (l42l , (I43K and (l44l give the far fields from the re- 
tardation, spatial-dispersion, and magnetic contributions, re- 
spectively. Despite the appearance of assorted phase factors 
at different stages of the derivation, to the extent that the re- 
fractive indices are real all net contributions have the same 
phase to within a plus or minus sign. The RD contributions in 
the two directions perpendicular to that of the polarization of 
the incident beam are equal, as expected by symmetry. This is 
not the case for SD and MG, since SD involves gradients and 
B is an axial vector. 
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As noted in the Introduction, the linear response cannot be 
calculated by factoring as done above. A full self-consistent 
Ewald-Oseen treatment is necessary. 



III. DISCUSSION 

Although Eqs. (l42l >. ( |43| >. and (l44l are complete, their gen- 
eral properties are not immediately evident. Hence we con- 
sider special cases. We also estimate conversion efficiency for 
fused silica, basing our calculations on several assumptions 
and the known LO properties of this material. 



A. Paraxial-ray approximation 

In the usual case of a highly collimated source beam of rel- 
atively small cross section, the emerging beam will also be 
initially relatively well localized but will diverge over a solid 
angle where the components essentially add in phase. Taking 
the diameter of the cross section W to be equal to at least a few 
wavelengths of the emerging beam, we make the paraxial-ray 
approximation, writing the observation direction for forward 
scattering as k = x6 x + y9 y + z, where the beam-divergence 
(observation) angles 6 X and 8 y are first-order quantities. With 
this representation the various projection operations are easily 
evaluated and we find 



[(/ — kk) ■ x](x ■ k) = xd x \ 

[(I-kk)-y](y-k) = ye y ; 
(I — kk) ■ z = —x6 x — yOy. 



(45) 



Considering also Eqs. d42li-(l44li. it is apparent that all contri- 
butions vanish in the forward direction and exhibit two-lobed 
patterns characteristic of gradient effects. Note that the z com- 
ponent also contributes on the same first-order scale when the 
viewer is off-axis. We shall use these equations in the follow- 
ing. 



B. Reduction to the isotropic case for large ki 

If the polarizable points are isotropic then Ci = Ct = C 
and Di = D t = D. If we assume further that Ki 3> 4mw 2 
then C w D, For clarity we write 



g(k,r, t) 



ttuj 2 NW 2 E 2 C 2 -(k 2 +k 2 )w 2 /s ihr-iiut 
8^ 6 



(46) 

since this is a common factor for all cases discussed in the rest 
of Sec.|m] Then Eqs. (I42t-(l44l reduce to 



E 



ff 

RD 



-2xn 2 x 



9(k,r,t) _ 
(n 2 - ni)' 



(47) 



E 



% = 4 [ xni 9 x - y{n 2 - m)0 y ] f k ^\ - (49) 

(n 2 — ni) 



where we have used the fact that k z differs from k only by 
terms of second order in 6. The net result is 



E 



ff 

Net 



2 [x (2ni + n 2 ) 6 X - y2(n 2 - n{)0 y ] 



g{k,r,t) 



{n 2 



m) 
(50) 

This limit applies to the SHG response of systems where the 
bond charge is strongly bound, for example organic materi- 
als and glasses. All three mechanisms generate polarization 
in the direction x of the applied field, and all have similar 
magnitudes, so none can be neglected. The fact that the RD 
term is important may not be at variance with the conclusion 
of Brudny et al.— which pertains to a configuration where the 
anharmonic contribution does not vanish completely. In par- 
ticular, the RD contribution here is exactly half that of SD and 
with opposite sign, so the net effect of the RD/SD combina- 
tion is to reduce the SD contribution by half. For n\ sa n 2 
the magnitude of the forward-scattered field intensity clearly 
benefits significantly from a long coherence length. 

Equation d49l shows that two y contributions are present, 
but to the extent that n\ w n 2 the overall term is small and 
can easily be overlooked since detection depends on intensity, 
not fields. This near-cancellation is a result of the sign of the 
z contribution in off-axis viewing. The cancellation of the 
y component is exact in Eq. d42l i. even in the general case 
where C\ ^ C t . A near-cancellation of the y component also 
occurs in the general case for Eq. (03]), although a second near- 
cancellation contributes if C\ ~ C t or D\ w D t . Thus if 
the y component is analyzed quantitatively, the more general 
equations must be used. We conclude that the z component is 
important in determining the properties of the emerging beam. 



C. Metals 

A second limit of the above is that corresponding to those 
metals for which the effective mass of the carriers is itself es- 
sentially isotropic. As a result of strong attenuation of op- 
tical signals, metals are usually measured in backscattering, 
where the results can be further complicated by surface recon- 



structions 



15.16.17 



Although much attention has been focused 
on these surface contributions j 17,18 ' 19 - 20 we consider here only 
signals originating in the bulk. In absorbing media the Green 
function retains its form, so the above development is still 
valid although the propagation vectors are now generally com- 
plex. With no restoring force K\ = 0, so C = AD. In the 
paraxial-ray approximation the RD contribution is unchanged, 
but the SD and MG terms are reduced by a factor of 4. For for- 
ward scattering the equations are 



E 



ff 

RD,m 



-2xn 2 9 x 



g(k,r,t) _ 
(n 2 - m) ' 



(51) 



E 



ff 

SD 



4X77,9 



g(k,r,t) _ 
(n 2 - m) 



(48) 



E 



■ff 

SD.m 



xn 2 t 



g(k,r,t) _ 
(n 2 - ni) ' 



(52) 
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E 



ff 

MG.m 



[xniO x - y(n 2 - ni)0 y 



g(k,r,t) _ 
(n 2 - nxY 



(53) 



(ni + n 2 J 



// 

TVet ■ 



-(n 2 - [z0 x - yflj,] 



g(k,r,t) 
(n 2 - ni) 



(54) 



The a; and y components now have equal amplitudes, but to 
the extent that n\ w n 2 the net result shows that the enhance- 
ment of the signal strength that results from the nearly singular 
denominator is cancelled. Again, all three contributing mech- 
anisms are important. Thus the assumption that the SHG con- 
tribution from metals arises entirely from spatial dispersion 
and magnetic effects is not quite correcti-t^^ 

SHG signals from surface reconstructions could be de- 
scribed in the above formalism by assigning suitable 
anisotropies to electrons in the surface region, although we 
do not do this here. 



D. Backscattering 

For backscattering the major difference is the reduction of 
the correlation length and corresponding reduction in the ra- 
diated field, since for negative k z the two terms of (2fc D — k z ) 
add instead of subtract. As we shall show in Sec. IIIIEl this 
effectively eliminates any possibility of observing SHG from 
the bulk of amorphous materials. The other effect is to reverse 
the sign of the result of the projection operation on z. When 
everything is taken into account, the paraxial-ray expressions 



for Ki ^> Amuj are 



E 



ff 

RD,b ~ 



2xn 2 8 x 



d(k,r,t) 
(ni + n 2 ) ' 



E 



ff 

SD,b 



-Axn 2 8 x x 



g(k,r,t) 
(ni + n 2 ) 



Emg b = 4[ot,i0 x + y{ni + n 2 )9 y ] 



g{k,r,t) 
(ni + n 2 ) 



(55) 



(56) 



(57) 



E f N {t, mb = K + n 2 )[x6 x + Wv] ? {k, *'*\ . (62) 

[ni + n 2 ) 



E. Power and conversion efficiency 

In many experiments what is determined is not the SHG 
intensity but the integrated SHG power. To obtain an order- 
of-magnitude estimate we consider the net x-polarized com- 
ponent for forward scattering with Ki ^> 4mw 2 and with 
Ci = C t = C and Dj = D t = D. The SHG intensity is 
given by 



Ish 



^l\E ff I 2 



(63) 



The SHG power is obtained by integrating this expression 
over a hemisphere of radius r. We are also interested in the 
conversion efficiency 77, which we define as 



n 



P S H 

(P "I 2 
y- 1 inc ) 



(64) 



where P- mc is the power of the incident beam. Assuming that 
the incident beam is collimated, the evaluation of its power in 
terms of the beam properties is straightforward, and we obtain 



P inc = r dx r dy C ^\E fe-^» w2 

J-00 J-00 o 71 " 



CTlf n, ,1 

16 



(65) 



That for the emerging beam is more complicated. The first 
issue concerns angular dependences. If the incident beam is 
reasonably well collimated and its diameter is equal to at least 
several SHG wavelengths, the SHG beam is also fairly well 
collimated. Then a small-term expansion in 8 is a good ap- 
proximation. To show this we consider 



E f Net,b = 2[g(2ni-na)fl z +2y(ni+n a )g v ] ff+^l • ( 58 > 

While both polarizations are present, the dominant contribu- 
tion in backscattering is that perpendicular to that of the driv- 
ing field. Although an x contribution is still generated, its 
strength is expected to be small compared to that polarized 
along y. 

The expressions for isotropic metals are 



-(kl+k 2 J/8W 



-(k 2 sin 2 6)/SW 



E 



ff 

RD.mb 



E 



ff 

SD,mb 



2xn 2 8 x 



-xn 2 ti x 



g(k,r,t) _ 
(m + n 2 ) ' 

g(k,r,t) 
(ni +n 2 )'' 



(59) 



(60) 



(66) 



Taking k = 2irn 2 /\sH, n 2 = 1.3, \sh = 400nm, and an in- 
cident beam width of 5 fim, we have k 2 /8W « 50. Hence the 
small-term approximation sinO ~ 8 is acceptable. This also 
provides justification for our use of the paraxial-ray approxi- 
mation in the previous sections. With these simplifications the 
area integral is straightforward and we find for x polarization 

cn 2 7rW^ 4 iV 2 £ 4 C 4 , ^ , 2 
Pru.x = — —-T-, — — (2ni + 712) 



8- 



16c 4 (n 2 - nif 

2tt 



; / d4> 8d8{8 2 cos 2 (, 
'0 Jo 



)e -k 2 e 2 w 2 /i (67) 



Tr 2 cq 2 N 2 E*C 2 
64n 2 (n 2 — n\) 2 



(2ni +n 2 f 
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Combining the above expressions we find the corresponding 
conversion efficiency to be: 



IV. CONCLUSIONS 



47: 2 q 2 N 2 C 4 (2n 1 +n 2 y 



(68) 



The efficiency decreases as the fourth power of the diameter 
of the incident beam. This is in contrast to the intensity, which 
decreases as 1/W 6 . 

From the definition of C; we have 



p = aEi oc = qAr = qCiEi oc , 



(69) 



where Ei oc is the field at the charge site and a is the linear 
polarizability. Then we can write C = a/q. We can connect 
a to the dielectric function t\ = n\ and bond density N of 
the material by the Clausius-Mossotti relation 



3 ei +2 



Then 



81(2ni + n 2 ) 2 
64n 2 cq 2 N 2 W 4 n 2 122(712 — n 



f-i + 2 



(70) 



(71) 



Using a driving wavelength A = 800 nm, dielectric functions 
of quartz of 2.112 and 2.161 at 800 and 400 nm, respec- 
tively, a bond density of 1.06 x 10 23 cm~ 3 f2i and a Gaus- 
sian beam of characteristic dimension W = 10/im we find 



T] X = 1.4 x 10" 



watt 



Thus 1 watt input power at 800 



nm is expected to generate about 3 SHG photons/sec. If W is 
reduced to 1/im, the output would increase to about 10 4 SHG 
photons/sec. These results appear to be consistent with exper- 
iment^ where few if any photons were seen emerging from 
the glass substrate. 



We have developed an anisotropic bond model (ABM) that 
describes SHG on the atomic scale, uses physically mean- 
ingful parameters, and includes all contributing mechanisms, 
thereby providing a more complete understanding of the 
physics of SHG than previously available. In disordered mate- 
rials the anharmonic restoring force acting on the bond charge 
does not contribute to the overall SHG signal, which instead 
arises from a combination of LO and beam-geometry effects 
and therefore provides limited new information about the ma- 
terial. For a Gaussian driving beam we obtain analytic expres- 
sions that give the phase, amplitude, and spatial distribution of 
the SHG radiation field for each of the remaining contributing 
mechanisms: retardation (RD), spatial dispersion (SD), and 
magnetic-field (MG) effects. All have the same order of mag- 
nitude, so any complete description must consider each. The 
expressions are reduced to simpler forms for both forward- 
and backscattering configurations in two isotropic limits, the 
first where the linear restoring force dominates, as in glasses, 
and the second where it is absent, as in metals. We esti- 
mate the conversion efficiency for forward scattering in fused 
quartz. Predictions appear to be in agreement with observa- 
tions, where available! 1 . Specific additional predictions allow 
critical tests of these results. 



With the basic physics established, we can now consider 
more complicated configurations, including nanospherical in- 
clusions in glass and the reported SHG enhancement with 
crossed-beam, crossed-polarization driving fields/ 1 " The re- 
sults presented here are also expected to be useful for ana- 
lyzing SHG data of liquids and biological materials. 
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